module gaussquad
	use ML
	implicit none
	integer, parameter	:: nGQ=100
	real				:: GQxw(nGQ,2)

	contains

	subroutine initgq()
		character(len=100)	:: nGQstring
		write (nGQstring,*) nGQ
		call execinML("GQxw=load('H:/mystuff/fortran/quadw"//trim(adjustl(nGQstring))//".dat')")
		call getML(GQxw,'GQxw')
	end subroutine
end module
    
module myfuncs
	use linds_int
	use lftds_int
	use lfdds_int

	use rnopt_int
	use rnnoa_int
	use rnunf_int
	use rnund_int
	use rnstt_int
	use rnset_int
	use rnfdf_int
	use rnun_int
	
	use tdf_int
	use tin_int
	use fin_int
	use anorin_int
	use anordf_int
	use eqtil_int
	use ordst_int
	use svrgn_int

	use uminf_int
	use bconf_int
	use erset_int
	use u4inf_int
	
	use qdval_int
	use bsnak_int
	use bscpp_int
	use csdec_int
	use pp1gd_int
	use bsitg_int
	use bs1gd_int
	
	use gamit_int
	use bsi0_int
	
	use ML
!	use imsl_libraries
	integer	time0, time1
	real, parameter		:: pi=3.14159265358979323846
	
!	private
	
!	public rann,invertpd,choleski,detpd,logdetpd,ransubsample,logmeanexp

    interface rann
        module procedure rann_m, rann_v, rann_s
	end interface

    interface setbounds
        module procedure setbounds_m,setbounds_v
	end interface

contains

	function gausscdfinv(x) result(val)
		implicit none
		real	:: x,val,a(1),y(1)
		a=x
		call vdcdfnorminv( 1, a, y )
		val=y(1)
	end function
	
	function gausscdf(x) result(val)
		implicit none
		real	:: x,val,a(1),y(1)
		a=x
		call vdcdfnorm( 1, a, y )
		val=y(1)
	end function
	
	subroutine setbounds_m(m,l,u)
		implicit none
		real	:: m(:,:),l,u
		where(m<l) 
			m=l
		elsewhere	
			where(m>u) m=u
		endwhere
	end subroutine
	
	subroutine setbounds_v(v,l,u)
		implicit none
		real	:: v(:),l,u
		where(v<l) 
			v=l
		elsewhere	
			where(v>u) v=u
		endwhere
	end subroutine

	subroutine inittime()
		call system_clock(time0)
	end subroutine
	
	subroutine printtime()
		character(len=1024) output
		
		call system_clock(time1)
		write (output,*) (time1-time0)/10000.0
		call VSprint(output)
	end subroutine
	
	elemental function boole(flag) result(val)
		logical, intent(in)	:: flag
		integer	:: val
		if(flag) then 
			val=1
		else
			val=0
		endif
	end function
	
	elemental function rboole(flag) result(val)
		logical, intent(in)	:: flag
		real	:: val
		if(flag) then 
			val=1
		else
			val=0
		endif
	end function

	function itostring(i) result(val)
		implicit none
		integer		:: i
		character(len=:),allocatable :: val
		character(len=100)	:: string
		write (string,*) i
		allocate(character(len=len(trim(adjustl(string))))::val)
		val=trim(adjustl(string))
	end function
		
    function getquadmin(x,y) result(val)
        implicit none
        real    :: x(3),y(3),val
        real    :: dy(3)
        
        dy(1)=y(2)-y(1)
        dy(2)=y(3)-y(1)
        dy(3)=y(3)-y(2)
        
        val=.5*(x(3)**2*dy(1)-x(2)**2*dy(2)+x(1)**2*dy(3))/(x(3)*dy(1)-x(2)*dy(2)+x(1)*dy(3))
        
    end function
    
    function getquadint(x,y,x0) result(val)
        implicit none
        real    :: x(3),y(3),x0,val
        
        val=((x0 - x(3))*((x0 - x(2))*(x(2) - x(3))*y(1) + (x0 - x(1))*(-x(1) + x(3))*y(2)) + (x0 - x(1))*(x0 - x(2))*(x(1) - x(2))*y(3))/((x(1) - x(2))*(x(1) - x(3))*(x(2) - x(3)))
        
    end function
    
    function getlinint(x,y,x0) result(val)
        implicit none
        real    :: x(2),y(2),x0,val
        
        val=((x0-x(2))*y(1)+(x(1)-x0)*y(2))/(x(1) - x(2))
        
	end function

	elemental function logit(x)
		real, intent (in)	:: x
		real				:: logit
		logit=exp(x)/(1+exp(x))
	end function

	function invertpd(A)
		real::A(:,:)
		real::invertpd(size(A,1),size(A,2))
		
		call linds(A,invertpd)
	end function
	
	function choleski(A) result(val)
		real::A(:,:)
		real::val(size(A,1),size(A,2))
		integer	:: j1,j2
		call lftds(A,val)
		do j1=1,size(A,1)
			do j2=j1+1,size(A,1)
				val(j1,j2)=0
			enddo
		enddo
	end function
	
	function diagonal(A) result(val)
		real::A(:,:),val(size(A,1))
		integer :: i
		do i=1,size(A,1)
			val(i)=A(i,i)
		enddo
	end function
	
	function detpd(A)
		real::A(:,:)
		real::detpd
		real::F(size(A,1),size(A,2)),det1,det2
		
		call lftds(A,F)
		call lfdds(F,det1,det2)
		detpd=det1*10**det2
	end function

	function logdetpd(A)
		real::A(:,:)
		real::logdetpd
		real::F(size(A,1),size(A,2)),det1,det2
		
		call lftds(A,F)
		call lfdds(F,det1,det2)
		logdetpd=log(det1)+log(10.0)*det2
	end function
	
	function rann_m(n,m)
		integer::n,m,i
		real::rann_m(n,m)
		
		do i=1,m
			call rnnoa(rann_m(:,i))
		enddo
	end function		
		
	function rann_v(n)
		integer::n
		real::rann_v(n)
	
		call rnnoa(rann_v)
	end function

	function rann_s()
		real::rann_s,dummy(1)
		call rnnoa(dummy)
		rann_s=dummy(1)
	end function

	function ranu_s() result(out)
		real::out
		out=rnunf()
	end function
	
	function ranu_v(n) result(out)
		integer::n
		real::out(n)
		call rnun(out)
	end function
	
	function ranui_v(n,k) result(v)
		integer		::n,k
		integer		::v(k)
		call rnund(n,v)
	end function
	
	function ranui_s(n) result(v)
		integer		::v,dummy(1)
		integer		::n
		call rnund(n,dummy)
		v=dummy(1)
	end function

	function logmeanexp(v) result(out)
		real		::v(:)
		real		::out,maxv
		maxv=maxval(v)
		out=maxv+log(sum(exp(v-maxv)))-log(real(size(v)))
	end function

	function logsumexp(v) result(out)
		real		::v(:)
		real		::out,maxv
		maxv=maxval(v)
		out=maxv+log(sum(exp(v-maxv)))
	end function

	function quantile(X,ps) result(qs)
		real		::X(:),ps(:),qs(size(ps)),xlo(size(ps)),xhi(size(ps))
		integer		nmiss

		call eqtil(X,size(ps),ps,qs,xlo,xhi,nmiss)
	end function
	
	function quantile_s(X,ps) result(qs)
		real		::X(:),ps,qs,xlo(1),xhi(1),psv(1),qsv(1)
		integer		nmiss

		psv(1)=ps
		call eqtil(X,1,psv,qsv,xlo,xhi,nmiss)
		qs=qsv(1)
    end function
    
	function stddev(X) result(std)
		real		::X(:),std
        std=sqrt(sum((X-sum(X)/size(X))**2)/size(X))
	end function
    
    function selectifc(X,cond) result(val)
		implicit none
		real	:: X(:,:)
		logical	:: cond(:)
		integer	:: i,c
		real, allocatable	:: val(:,:)
		if(size(X,2).ne.size(cond)) then
			call VSprint("array size don't match in selectifc")
			stop
		endif
		allocate(val(size(X,1),count(cond)))
		c=1
		do i=1,size(X,2)
			if(cond(i)) then 
				val(:,c)=X(:,i)
				c=c+1
			endif
		enddo
	end function

	function selectifr(X,cond) result(val)
		implicit none
		real	:: X(:,:)
		logical	:: cond(:)
		integer	:: i,c
		real, allocatable	:: val(:,:)
		if(size(X,1).ne.size(cond)) then
			call VSprint("array size don't match in selectifr")
			stop
		endif
		allocate(val(count(cond),size(X,2)))
		c=1
		do i=1,size(X,1)
			if(cond(i)) then 
				val(c,:)=X(i,:)
				c=c+1
			endif
		enddo
	end function
		
end module	
		
		
